日本の新型コロナウィルスに関する情報が抱える大きな問題は、ソースごとに値が微妙にことなっていることです。原因は明白で、総元締めであるはずの厚生労働省が各自治からの患者データを集計せずに各自治体のプレスリリースやウェブサイトの数値を積上げたものを発表していることにあります。
本来であれば医療機関->市町村->都道府県->厚生労働省の順で患者のデータを全て吸い上げて集計すべきなのです。

Import

厚生労働省

厚生労働省のオープンデータのページから取得した「陽性確定者数」のデータ

df_mhlw <- "https://www.mhlw.go.jp/content/pcr_positive_daily.csv" %>% 
  readr::read_csv() %>% 
  dplyr::mutate(date = lubridate::as_date(`日付`),
                n = as.integer(`PCR 検査陽性者数(単日)`),
                cum = cumsum(n)) %>% 
  dplyr::select(date, n, cum)

df_mhlw

 

Covid19Japan

Covid19 Japanが独自に収集している陽性確定者単位のデータ。ソースとデータは全てGitHubにて公開されているが、データはJSON形式である点に注意。発表後に修正されたレコード(インスタンス)は削除されれずにステータスなどが変更されているだけなので、「レコード数 \(\neq\) 累計陽性確定者」である点に注意。

df %>% 
  dplyr::group_by(patientId) %>% 
  dplyr::summarise(n = n())
df %>% 
  dplyr::filter(patientId == "-1")
df %>% 
  dplyr::filter(patientStatus == "Deceased")

JAG Japan

ジャッグジャパン株式会社が独自に収集している陽性確定者単位のデータ。データはCSV形式にて公開されているが、GIS処理用のデータが含まれている点に注意。「レコード数 \(=\) 累計陽性確定者数」になっている。

df_jag <- "https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv" %>% 
  readr::read_csv()
## Warning: Missing column names filled in: 'X52' [52], 'X53' [53], 'X54' [54]
## Warning: 610 parsing failures.
##   row              col expected              actual                                                               file
##  4484 市町村内症例番号 a double 神戸59/西宮36       'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
##  4485 市町村内症例番号 a double 神戸60/西宮37       'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
##  6054 市町村内症例番号 a double 川口40/さいたま50   'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 12486 市町村内症例番号 a double 神戸237/西宮63      'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 14454 市町村内症例番号 a double (271)のちに番号なし 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## ..... ................ ........ ................... ..................................................................
## See problems(...) for more details.
df_jag

 

Area Data

都道府県と地方区分に関するデータをGistで公開していますので、使いたい方はご自由に。

 

Data Wrangling

ソースによる違い

収集組織 累計感染者数 累計死者数 日付
厚生労働省 96948 2020-10-26
Covid19 Japan 98146 1726 2020-10-27
JAG Japan 93661 2020-10-26
NHK(参考) 98262 1733 2020-10-27

https://www3.nhk.or.jp/news/special/coronavirus/data-all/

Summarize

オリジナルのデータがどのようになっているかskimrパッケージを用いてサマライズしてみます。元がJSON形式なので、読み込んだ直後は殆どの変量(フィーチャー)が文字型になっているのと欠損値が多いことが分かります。

df %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 99726
Number of columns 23
_______________________
Column type frequency:
character 19
logical 3
numeric 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
patientId 0 1.00 1 8 0 98147 0
dateAnnounced 0 1.00 10 10 0 273 0
gender 12662 0.87 1 1 0 2 0
detectedPrefecture 0 1.00 3 15 0 49 0
patientStatus 95831 0.04 8 23 0 8 0
notes 51724 0.48 1 270 0 45257 1
mhlwPatientNumber 99277 0.00 1 11 0 434 0
prefecturePatientNumber 10702 0.89 5 20 0 89015 0
prefectureSourceURL 68470 0.31 5 224 0 3433 0
residence 20455 0.79 1 38 0 1421 0
sourceURL 637 0.99 1 239 0 7736 0
relatedPatients 89684 0.10 2 259 0 6167 0
knownCluster 97269 0.02 3 88 0 228 0
detectedCityTown 74545 0.25 2 22 0 663 0
cityPrefectureNumber 74817 0.25 1 34 0 24900 2
citySourceURL 87930 0.12 9 317 0 3623 0
deceasedDate 97990 0.02 10 10 0 223 0
deceasedReportedDate 98609 0.01 10 62 0 183 0
deathSourceURL 98706 0.01 14 123 0 624 0

Variable type: logical

skim_variable n_missing complete_rate mean count
confirmedPatient 0 1 0.98 TRU: 98146, FAL: 1580
charterFlightPassenger 99712 0 1.00 TRU: 14
cruisePassengerDisembarked 99715 0 1.00 TRU: 11

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
ageBracket 0 1 32.98 23.44 -1 20 30 50 100 ▃▇▅▂▁

 

Transform

各変量(フィーチャー)を適切な形式に変換し、地域区分でも分析できるように都道府県データと結合します。

x <- df %>% 
  dplyr::mutate(dateAnnounced = lubridate::as_date(dateAnnounced),
                ageBracket = forcats::as_factor(ageBracket),
                gender = forcats::as_factor(gender),
                patientStatus = forcats::as_factor(patientStatus),
                residence = forcats::as_factor(residence),
                cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE)) %>% 
  dplyr::select(dateAnnounced, ageBracket, gender, detectedPrefecture,
                patientStatus, knownCluster, cluster) %>% 
  dplyr::left_join(prefs, by = c("detectedPrefecture" = "pref"))

x
x %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 99726
Number of columns 13
_______________________
Column type frequency:
character 2
Date 1
factor 9
logical 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
detectedPrefecture 0 1.00 3 15 0 49 0
knownCluster 97269 0.02 3 88 0 228 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dateAnnounced 0 1 2020-01-15 2020-10-27 2020-08-10 273

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
ageBracket 0 1.00 FALSE 13 20: 23833, 30: 14960, -1: 12773, 40: 12330
gender 12662 0.87 FALSE 2 M: 48990, F: 38074
patientStatus 95831 0.04 FALSE 8 Dec: 1726, Hos: 1266, Hom: 316, Dis: 283
pcode 1159 0.99 FALSE 47 13: 30790, 27: 12272, 14: 8559, 23: 6031
都道府県 1159 0.99 FALSE 47 東京都: 30790, 大阪府: 12272, 神奈川: 8559, 愛知県: 6031
八地方区分 1159 0.99 FALSE 8 関東地: 52254, 近畿地: 19525, 九州地: 10905, 中部地: 9643
広域圏 7448 0.93 FALSE 8 首都圏: 52471, 近畿圏: 18965, 中部圏: 8228, 九州圏: 7592
通俗的区分 1159 0.99 FALSE 11 関東: 52254, 関西: 18965, 東海: 7882, 九州: 7592
fct_pref 1159 0.99 FALSE 47 Tok: 30790, Osa: 12272, Kan: 8559, Aic: 6031

Variable type: logical

skim_variable n_missing complete_rate mean count
cluster 0 1 0.02 FAL: 97269, TRU: 2457

 
文字型を因子型に変換するだけでも、変量(フィーチャー)内の水準の比率が大まかにつかめるようになります。
例えば、年代別の陽性判定者数は20代が最も多く、続いて、30代、40代と働き盛りの世代に多いことが分かります。また、都道府県別では、東京、大阪、神奈川、愛知と成っていますが、地方区分別では、以外にも関東、大阪の次に九州がきていることが分かります。

 

Visualizing

厚生労働省

df_mhlw %>% 
  ggplot2::ggplot(ggplot2::aes(x = date)) + 
    ggplot2::
                  ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)

全体傾向

sec_scale <- 50

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n),
                ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)) %>%
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
                      fill = "dark green", alpha = 0.5, width = 1.0) + 
    ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
                       colour = "dark green") +
    ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark red") + 
    ggplot2::scale_y_continuous(
      name = "陽性確定数(日別)・7日間移動平均(濃赤折線)",
      sec.axis = ggplot2::sec_axis(~ . * sec_scale,
                                    name = "累積陽性確定者数(折線)")
    ) +
    ggplot2::labs(title = paste0("Covid19, Japan(@", Sys.time(), ")"),
                  subtitle = "", x = "日付")
## Warning: Removed 6 row(s) containing missing values (geom_path).

plotly::ggplotly()

地方区分別

都道府県別では区分が多すぎるので八地方区分別の積上げ棒グラフを描いてみます。サマリで見たように九州地方が以外にも多いことが読み取れます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

折線グラフで表示してみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::geom_point(ggplot2::aes(y = n), size = 1) + 
    ggplot2::scale_colour_brewer(palette = "Set1")

 
さらに分かりやすくするために関東、中部、近畿、九州の四地区を除く他の地区をその他としてまとめてみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::scale_colour_brewer(palette = "Set1") +
    ggplot2::facet_wrap(~ `八地方区分`, ncol = 2)

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `八地方区分`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "八地方区分", values_to = "n") %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_relevel(`八地方区分`,
                                              "北海道地方", "東北地方",
                                              "関東地方", "中部地方",
                                              "近畿地方", "中国地方",
                                              "四国地方", "九州地方")) %>%
  dplyr::group_by(`八地方区分`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `八地方区分`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19(@", Sys.time(), ")"),
                  y = "累計人数")

 
九州が8月頃から急増しているのは、県別に見ると福岡と沖縄での急増が原因と分かります。

x %>% 
  dplyr::filter(`八地方区分` == "九州地方") %>% 
  dplyr::group_by(dateAnnounced, `都道府県`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `都道府県`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "都道府県", values_to = "n") %>% 
  dplyr::mutate(`都道府県` = forcats::fct_relevel(`都道府県`,
                                              "福岡県", "佐賀県", "長崎県",
                                              "熊本県", "大分県", "宮崎県",
                                              "鹿児島県", "沖縄県")) %>%
  dplyr::group_by(`都道府県`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `都道府県`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19, 九州地方(@", Sys.time(), ")"),
                  y = "累計人数")

 

クラスタ比率

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(`八地方区分`, cluster) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::pivot_wider(names_from = cluster, values_from = n) %>% 
  dplyr::mutate(ratio = (`TRUE` / (`TRUE` + `FALSE`) * 100) %>% round(1)) %>% 
  dplyr::rename(`地方` = `八地方区分`,
                `非クラスタ感染者` = `FALSE`, `クラスタ感染者` = `TRUE`,
                `クラスタ比率[%]` = ratio)

Modeling

時系列(TS)分析

日毎の陽性判定者数を時系列データ形式に変換します。判定者数がゼロの日は報告されていないので、最初の陽性判定者が出た日から日日次のカレンダーデータを作成して結合してます。なお、時系列データ形式の周期については日次データなので7日間を周期として設定しておきます(考え方によっては約1ヶ月ごとで4週間=28日を周期にするのもありかと思います)。

tmp <- x %>%
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% print() %>% 
  dplyr::rename(date = dateAnnounced) %>% 
  tidyr::complete(date = seq.Date(min(date), max(date), by = "day"),
                  fill = list(n = 0L)) %>% print() %>% 
  # tidyr::replace_na(list(n = 0L)) %>% print() %>% 
  dplyr::select(-date)
## # A tibble: 270 x 2
##    dateAnnounced     n
##    <date>        <int>
##  1 2020-01-15        1
##  2 2020-01-24        1
##  3 2020-01-25        1
##  4 2020-01-26        1
##  5 2020-01-28        3
##  6 2020-01-29        1
##  7 2020-01-30        4
##  8 2020-01-31        1
##  9 2020-02-01        2
## 10 2020-02-04        1
## # … with 260 more rows
## # A tibble: 287 x 2
##    date           n
##    <date>     <int>
##  1 2020-01-15     1
##  2 2020-01-16     0
##  3 2020-01-17     0
##  4 2020-01-18     0
##  5 2020-01-19     0
##  6 2020-01-20     0
##  7 2020-01-21     0
##  8 2020-01-22     0
##  9 2020-01-23     0
## 10 2020-01-24     1
## # … with 277 more rows
ts_x <- tmp %>%
  ts(frequency = 1)

tsw_x <- tmp %>%
  ts(frequency = 7)

tsm_x <- tmp %>%
  ts(frequency = 28)
# ts_x <- x %>% 
#   dplyr::filter(!is.na(fct_pref)) %>% 
#   dplyr::group_by(dateAnnounced) %>% 
#   dplyr::summarise(n = n()) %>% 
#   dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
#   dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
#   dplyr::select(-date) %>%
#   ts(frequency = 1)
#   
# tsw_x <- x %>% 
#   dplyr::filter(!is.na(fct_pref)) %>% 
#   dplyr::group_by(dateAnnounced) %>% 
#   dplyr::summarise(n = n()) %>% 
#   dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
#   dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
#   dplyr::select(-date) %>%
#   ts(frequency = 7)
# 
# tsm_x <- x %>% 
#   dplyr::filter(!is.na(fct_pref)) %>% 
#   dplyr::group_by(dateAnnounced) %>% 
#   dplyr::summarise(n = n()) %>% 
#   dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
#   dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
#   dplyr::select(-date) %>%
#   ts(frequency = 28)

時系列データに変換したものをプロットすると可視化の項でプロットした棒グラフと同じ形のグラフになることが分かります。

ts_x %>% 
  plot()

tsw_x %>% 
  plot()

tsm_x %>% 
  plot()

上記からトレンド(長期的傾向)を除いたグラフ。デフォルト指定なのでlag = 1

ts_x %>% 
  base::diff() %>% 
  plot()

tsw_x %>% 
  base::diff() %>% 
  plot()

tsm_x %>% 
  base::diff() %>% 
  plot()

トレンド、季節変動(周期変動)、非周期変動に分解した場合。frequency = 1では分解できない点に注意。

tsw_x %>% 
  stats::decompose() %>% 
  plot()

tsm_x %>% 
  stats::decompose() %>% 
  plot()

トレンドだけを抜き出してみる。

tsw_x %>% 
  stats::decompose() %>% 
  .$x %>% 
  plot(ylim = c(0, 1500))

par(new = TRUE)

tsw_x %>% 
  stats::decompose() %>% 
  .$trend %>% 
  plot(ylim = c(0, 1500), col = "dark green", lwd = 3)

par(new = FALSE)
tsm_x %>% 
  stats::decompose() %>% 
  .$x %>% 
  plot(ylim = c(0, 1500))

par(new = TRUE)

tsm_x %>% 
  stats::decompose() %>% 
  .$trend %>% 
  plot(ylim = c(0, 1500), col = "blue", lwd = 3)

par(new = FALSE)

# tsm_x %>% 
#   stats::decompose() %>% 
#   .$trend %>% 
#   plot()

スペクトル密度の推定。何を意味するのかよく分からない。

ts_x %>% 
  spec.pgram()

tsw_x %>% 
  spec.pgram()

tsm_x %>% 
  spec.pgram()

自己回帰によるスペクトラム解析

spectrum(ts_x, method="ar")

spectrum(tsw_x, method="ar")

spectrum(tsm_x, method="ar")

Infer

ARIMA Model

ARIMA(Auto Regressive Integrated Moving Average, 自己回帰和分移動平均)モデルによる予測を行ってみます。予測に必要なパラメータはステップワイズにより自動的に最適なものが選択されます。

ts_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsw_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsm_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()